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Abstract 

We discuss tiow semidefinite programming can be used to determine the second-order density 
matrix directly tlirough a variational optimization. We show how the problem of characterizing 
a physical or A^-representable density matrix leads to matrix-positivity constraints on the density 
matrix. We then formulate this in a standard semidefinite programming form, after which two 
interior point methods are discussed to solve the SDP. As an example we show the results of an 
application of the method on the isoelectronic series of Beryllium. 
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I. INTRODUCTION 



The idea of a variational determination of the ground-state energy for a non-relativistic 
many-body problem based on the second-order density matrix (2DM) has a long history [!}- 
[3j and several highly appealing features. The energy of a system is a known linear functional 
of the 2DM. A^-particle wave functions never need to be manipulated since the energy is 
minimized directly in terms of the 2DM. However, the minimization is constrained because 
the variational search should be done exclusively with 2DMs that can be derived from an 
A^-particle wave function (or an ensemble of A^-particle wave functions). Such a 2DM is 
called iV-representable, and the complexity of the many-body problem is in fact shifted 
to the characterization of this set of A^-representable 2DMs. The complete (necessary and 
sufficient) set of conditions for A^-representability of a 2DM is not known in a constructive 
form, but it is clear that the energy from a minimization constrained by a set of necessary 
A^-represent ability conditions is a strict lower bound to the exact energy. Therefore this 
approach is highly complementary to the usual variational procedure based on a wave- 
function ansatz, which produces upper bounds. In addition the method is in principle 
exact, in the sense that as increasingly accurate set of A^-representability conditions are 
imposed in the minimization, the resulting energy converges to the exact one. 

These are fascinating ideas for any true-blooded many-body theorist, as it comes close 
to the "ultimate reduction" of an interacting many-particle problem to solving a sequence 
of two-particle problems. In practice, however, implementing the method turns out to be 
very difficult and it is only in the last decade that serious attempts have been undertaken 
to turn the idea into a practical calculational scheme. The efforts by Mazziotti et al. [1H6] 
and Nakata et al. [TJ [8] are particularly notable. The main difficulty is of a technical na- 
ture: stringent A^-representability conditions require the positive semidefiniteness of matrix 
functionals of the 2DM, which turns the variational problem into a so-called semidefinite 
program (SDP). Even applying the simplest "two-index" conditions, a direct energy mini- 
mization using Newton- Raphson methods requires a matrix operation scaling as M^^ (where 
M is the number of single-particle states) in each Newton-Raphson step. This can be cir- 
cumvented in various ways, so that only matrix operations scaling as are needed. While 
these are nominally methods, the number of iterations required to reach convergence 
is very high and seems to rise with system size; in practice present implementations are 
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probably about 100-1000 times slower than comparable methods. Still, one has the feeling 
that there is potential to turn it into a genuine method, and it is of interest to inves- 
tigate the properties of SDP applied to various systems. In this proceeding we will first 
discuss the problem of density matrix optimization and A^-representability, and how it can 
be formulated as an SDP. We will describe two different algorithms we used to solve the 
problem and as an example we will discuss the results of application of the method to the 
isoelectronic series of Be [9] . 



II. VARIATIONAL DENSITY MATRIX DETERMINATION 

We use second quantized notation where aj^ (aa) creates (annihilates) a fermion in a 
single-particle (sp) state a). When there are only two-body interactions present in a physical 
system, the Hamiltonian of that system can be written as: 

The expectation value of the energy in an arbitrary A^-particle state |\E'^) can be expressed 
using only the second-order density matrix (2DM), 

E{r) = Tr ri7(2) = \ J2 r«/5;7^<;.. , (2) 

a/375 

with the 2DM defined as: 

r«/3;75 = (^^|a],a^a5a^|^^) , (3) 
and the reduced two-particle Hamiltonian, 

Now the idea of variational density matrix optimization is to determine the ground-state 
energy and other two-body properties by minimizing the energy (|2| using the 2DM as a 
variable. This is a much compacter object than the wavefunction because, no matter how 
many particles are involved, you always stay in two-particle (tp) space. The problem is 
that there is no straightforward way to know whether an arbitrary matrix in tp-space F is 
derivable from a physical A^-particle fermionic wavefunction as in Eq. (|3]). This is called 
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the A^-representability problem. Some obvious necessary A^-representability constraints are 
apparent from the definition ([3]): 

T^T^^. (5) 

^al3]-yS = — r^a;7(5 = ^ r^^^^^ = T pa-S-y , (6) 
^al3]-yS = r^(5;a^ , (7) 

but it turns out that there are many more non-trivial constraints needed to ensure that a 
2DM is physical. 



A. A^-representability 

The necessary and sufficient conditions for A^-representability are known. A tp-matrix is 
iV-representable if and only if, for every two-body Hamiltonian H^, the following inequality 
is satisfied: 

Tr Hi^^T > Eo{H,) . (8) 

This is of course not a constraint that can be used in practice, as you need to know the 
ground-state energy of every imaginable two-body Hamiltonian. Therefore one resorts to 
certain classes of Hamiltonians for which a lower bound to the ground-state energy is known 
[21 El Hn] • A Hamiltonian class that is used as necessary constraint is 

{^^\B^B\^^) > , (9) 

which leads to positivity conditions of linear matrix maps of the 2DM. There are three 
possible forms of the operator i?^ if we want (joj) to be expressable as a function of the 2DM, 
which gives rise to three conditions on the density matrix: 
a. B"^ = ^apPapCi^a'^^ Icads to the trivial X-condition: 

Y,Pap{'^''\aia\asa^\^'')p.,5>() ^ X(r) = T ^ , (10) 

which just demands that the 2DM remains positive semidefinite. 
ai3 Qai^aaap leads to the Q-condition: 

5^ga/3(^'^|a«a^4at|^^)g^5>0 ^ Q(r) ^ , (11) 

in which the Q can be written as a function of F using the anticommutation relations. 
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c. = J2ai3 QaiiCi^a'^p which leads to the ^-condition: 



^ ga,p{m''\aiapa\a^\^>'')g^s > -> ^^(r) h , (12) 

a/375 

Another Hamiltonian class for which a lower bound to the ground-state energy is known, 
gives rise to the so-called three-index conditions: 

>0 . (13) 



Two commonly used three- index conditions can be derived from Eq. (13): 

d. B"^ = J2ai3-f ^a/37^a^/3"7 ^^ads to the Ti-couditiou: 

J2 tlp^{^^\aiala\ac;a,as + a^a/ja^a^ajaj |^'^)4c > ^ TiiT) >z . (14) 

e. i?"!" = Y.al3y ^a/37'^a'^li'^7 ^^^'^^ tO the T^-COuditioU 

tlp^i^^laiala^ala.as + a^a^a^^aja^ |^^)4c > ^ T2{T) h . (15) 

In short, the optimization problem that we have to solve can be summerized as: 

mm Tr TH^^^ , (16) 

under the condition that 

ivr^^. (17) 

£(r)^o \/ce{i,Q,g,ruT2} . (is) 

III. REPRESENTATION AS AN SDP 

The variational method described in the previous section can be formulated as a semidef- 
inite program. A general 2DM, describing an A^-particle system can be expanded in an 
arbitrary orthogonal basis of traceless matrix space {/'} as 

with M the dimension of single-particle (sp) space, and the unit matrix on tp space defined 



as 



(^tp)a/3;76 ~ ^a-y^fSS " SaS^/Sj ■ (20) 
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The energy of the system can be written as a function of the 7's as 

Tr TH^^) = ^S^'^^ Tr H^'^ + V 7,Tr H^'^ f . (21) 
M(M- 1) i 

Because the necessary A^-representabihty conditions can be written as hnear homogeneous 
matrix maps of F, we can also write them as a function of the 7's: 

^ = M{L~-!f + ^ - ° • (2') 



If we now define the block matrices: 

N{N-l) 



u 



M(M- 1) 



0£,(ltp) and u^ = Qc,{r) , (23) 



then we can formulate the variational density matrix optimization problem as a standard 
dual form of a semidefinite program [TT] : 

min 7^/i on condition that Z = +y ^iU^ y , (24) 
7 



in which /i* = Tr H^"^^ f\ The primal problem corresponding to (24) optimizes the matrix- 
variable X, the problem is defined as: 

max - Tr Xm° on condition that Tr Xu' = h' and X >- . (25) 

The primal-dual gap r/ is defined as the difference between the primal and the dual cost 
function for a certain primal-dual point (X, Z\. 

T^ = ^'^h + Tr u^X = 7iTr Xu' + Tr Xu^ = Tr > , (26) 

i 

because X and Z are positive semidefinite matrices. We can see that the smallest value 
of 7] will be reached when both the primal and the dual problem are optimal. It can be 
proved [H] that if the primal and the dual problem are both feasible, then at their solution 
the primal-dual gap will vanish. This means that the primal-dual gap can be used as a 
convergence criterion for the algorithm, what is more, at any point during the optimization, 
the error on the current value is limited from above by the primal-dual gap. 



IV. ALGORITHMS USED TO SOLVE SDP 



There is a vast literature available with different methods to solve SDP's. It is generally 
excepted that the interior point methods have the best computational performance. In the 
course of our research we have implemented numerous algorithms, here we will discuss two 
of them, a dual-only potential reduction method, and a primal-dual path following method. 



A. Dual potential reduction method 



This method will only solve the dual problem (24), so we have no access to the primal- 
dual gap as a convergence criterion. The idea is to minimize the following potential function 
over 7: 

$(7, t) =j'^h-t In det Z(7) . (27) 

Starting from a feasible Z, the logarithmic potential will make sure that Z remains positive 
definite. The potential is minimized for decreasing scaling factor t, using the solution of 
the previous t as a starting point. When t — i- the solution will lie on the edge of the 



feasible area, and this is exactly the solution of the original semidefinite program (24). For 



every value of t the optimization problem of (27) will be solved using Newton method. This 



means that we have to solve the following linear system of equations to find the step ^7 until 
convergence is reached: 

t J2 (Tr u"Z-\^Z-^) {6-f)j = tTiZ-W - h' . (28) 
j 

The dimension of the system of equations is M'^, which means that if we would solve it by 
direct inversion the method would scale as M^^. We can construct an efficient matrix- vector 
product and use the conjugate gradient method to solve this much faster. Unfortunately in 
the limit of small t the number of iterations needed to converge increases because the system 
becomes ill conditioned. 



B. Primal-dual path following method 

In this method we will solve the primal and the dual problem at the same time. To 
explain how it works we first have to define the central path, which is the set of primal-dual 
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points for which 



XZ = ^ls,p , (29) 
n 



with n the total dimension of the X and Z matrices and Igup the direct sum of the unity 
matrices on the different constraint spaces: 

Isup = Ifc • (30) 
k 

In the path following algorithm we will try to follow the central path, reducing the primal- 
dual gap along the way. Imagine that we have a primal-dual point {X, Z) on the central 
path with primal-dual gap r]. We want to answer the question: what is the primal-dual point 
on the central path with primal-dual gap scaled down with a factor z/, rephrased, what are 
the (AX, AZ) that solve: 

(X + AX)(Z + AZ) = ^l,„p. (31) 

There are several ways to symmetrize these equations, using the method proposed by Nes- 
terov and Todd (see |12j), we obtain two equivalent equations, which we call the primal and 
the dual equation: 

(P) : AZ + D-^AXD-^ = '^X-^ - Z , (32) 



with 



(D) : AX + D AZ D = —Z'^ - X , (33) 

n 



D{X, Z) = Z--^ (z^^xz'^^ ' Z"^ , (34) 
and under the condition that: 

Tr AXu' = and AZ = ^(57)^^* . (35) 

i 

We can see that now, we will have to solve two systems of linear equations: 

J2 (Tr D c> D d) Szj = Tr (^X"^ -Z^d , 

j 

J2 (Tr D-\W-\') (57), = Tr (^^"' - u\ 



n 



in which we define {d} as the basis of the orthogonal complement of the space spanned by 
the u*'s, and AZ = ^j^Zjd. We will again be able to construct an efficient matrix-vector 
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product and use the conjugate gradient method. First we solve the dual system, which has 
the smallest dimension, and feed the solution into the primal system, which then needs very 
few iterations to converge. Once again we run into the same problem that the solution of 
the dual system needs more and more iterations as it gets closer to the solution, because of 
the increasing condition number of the system. 

V. APPLICATION TO BERYLLIUM ISOELECTRONIC SERIES 

We have applied this method in the study of electronic structure calculations of diatomic 
molecules and atomic systems [9j. We will discuss in short the result of the method 

applied to the isoelectronic series of Beryllium, using the P, Q and G conditions. This is 
an interesting test case for the inclusion of static electron correlation, because there is a 
near degeneracy between the 2s and 2p energy levels when the central charge becomes large, 
thus creating a multireference ground state. Because of the amount of symmetry present in 
atomic systems, the density matrix and the constraint matrices all become block diagonal 
(see [9]), enabling us to use quite large basis sets (cc-pVDZ, cc-pVTZ and cc-pVQZ [TU]). 
In Fig. [T] the SDP correlation energy, defined as the difference between the Hartree-Fock 
and the SDP energy, is shown as a function of central charge Z for the different basis sets, 
compared to estimates for non-relativistic energies based on experimental data [171 118] , and 
to the results of coupled cluster (CCSD) calculations. The experimental curve is linear in Z, 
as a direct consequence of the near-degeneracy of the ground state [18j. The SDP correlation 
energy does not follow this trend: it goes linear in the beginning, but becomes concave in 
the cc-pVDZ and cc-pVTZ basis, or convex in the cc-pVQZ basis. This failure, however, 
is not related to the SDP method as the trend is the same in the CCSD calculations. It 
simply reflects the fact that the incipient degeneracy is not well described in these basis sets. 
This can also be seen by calculating the hydrogen spectrum (corresponding to the Z ^ oo 
situation, when the electron-electron interaction can be neglected) in the basis sets: the 
2s and 2p energies are not degenerate, but differ by 5.8 mhartree (cc-pVDZ), 2.0 mhartree 
(cc-pVTZ) and -2.3 mhartree (cc-pVQZ). We also performed calculations in the cc-pVDZ 
basis after rescaling (r — )■ ar) it in such a way that the hydrogenic 2s-2p degeneracy is exact. 
In this basis the SDP correlation energy (also shown in Fig. [T]) indeed has the correct linear 
behavior. It is clear from the above discussion that SDP is indeed capable of providing 
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FIG. 1. The SDP correlation energy (£^hf — -E^sdp) for the Be series in all three basis sets, and in 
a rescaled basis set that exhibits hydrogen-like behavior (degeneracy between the 2s and 2p level) . 
For comparison, the CCSD and experimental values are also shown. 

accurate correlation energies in the presence of near degeneracies, when other many-body 
techniques (like density functional theory or MP2) can fail. 
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